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Abstract. This contribution summarizes the recent work carried out to analyze the behavior 
of the hyperbolic sector of the Fully Constrained Formulation (FCF) derived in Bonazzola et al. 
2004. The numerical experiments presented here allows one to be confident in the performances 
of the upgraded version of CoCoNuT's code by replacing the Conformally Flat Condition (CFC) 
approximation of the Einstein equations by the FCF. 



1. Introduction 

Most numerical codes used to solve the Einstein equations in order to obtain stationary 
or dynamical spacetimes generated by compact astrophysical objects are based on the 3+1 
formalism (see, e.g., [T]). The Fully Constrained Formalism (FCF), recently proposed in [21[31ll], 
is a constrained evolution formulation of the full Einstein equations. It is a natural generalization 
of the Conformally Flat Condition (CFC) [51 E], which is an approximation of the Einstein 
equations. FCF extends the CFC approximation and includes a hyperbolic system governing 
the gravitational radiation. Reader interested in details and motivation about FCF can address 
to the reference [2]. We present numerical simulations of this system in several cases. 



2. Formalism 

Given an asymptotically flat spacetime {AA^g^y) we consider a 3 + 1 splitting by spacelike 
hypersurfaces Sj, with being the timelike unit normal to Tit- Latin (greek) indices go from 
1 (0) to 3. '^^u = g^p + n^Ui, denotes the 3-metric on and K^y = —^Cn^y^u the extrinsic 
curvature. With the lapse function N and the shift vector the metric g^p is expressed in 
coordinates {x^^) as g^y dx^ dx" = —N"^ dt^ + ^ij{dx^ + /3* dt){dx^ + fi^ dt). As in [2], we introduce 
a flat metric /jj, which satisfies dtfij = and fij ~ 7ij at spatial infinity. We define 7 := det7jj 
and / := det/jj. We introduce the following conformal decomposition: 7^ = ip^^ij. We define 
f^ij .— ^ij _ ji] _ fjij^g gauge in is maximal slicing, K = 0, and the so-called generalized Dirac 
gauge, T^kl^^ = 0, where Pfc is the Levi-Civita connection associated with fij. The Einstein 
equations become a coupled elliptic-hyperbolic system: the elliptic sector acts on the variables 
"0, A^, and and the hyperbolic sector acts on /i*-' [2j . If /i*-' = is imposed, CFC is recovered. 
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Figure 1. Radial profile of h^'^ at t = 6. Right figure corresponds to the equator and left figure 
corresponds to the pole. Continuous lines corresponds to = 50 and ne = 100 (radial and 
angular number of grid points), = 100 and ng = 20, and = 200 and ng = 40, in non-linear 
evolutions. Dashed line corresponds to n,. = 100 and ng = 20, in a linear evolution. 



We introduce the conformal decomposition A^^ = ip^^ (^K^^ — ^K^^^^, where K = j^^Kij. 
This decomposition is different from the one introduced in [2], but it is motivated by the local 
uniqueness properties of the elliptic equations shown in [3]. We rewrite the equation for h^^ as 
a first order evolution system for the tensors {h^^ , A^^ , w^^ ) , where w^^ := I'fcT*-' [7J. 

3. Numerical simulations 

We evolve numerically the previous evolution system for (h^^ , A^^ , w^^ ) . During the numerical 
simulations, we consider N, f3'^, ip and the energy-momentum tensor as sources of the system. 
We perform the evolution of matter with the CoCoNuT code [S]. Some basic elements of this 
code are: i) fourth order finite differences scheme for the spatial derivatives and fourth order 
Runge-Kutta methods for the time derivative; ii) axisymmetry and symmetry with respect to 
the equatorial plane; iii) spherical orthonormal coordinates, {r,9,cp); iv) a Sommerfeld condition 
at the outer boundary; and v) Kreiss-Oliger dissipative term in order to avoid the numerical 
noise of high frequency which appears during long-term simulations. 

3.1. Teukolsky waves 

The first test is the evolution of a combination of ingoing and outgoing even-parity axisymmetric 
Teukolsky waves [9j, in order to construct regular initial data at the center r = 0. The amplitude 
chosen is 10~^. These data are a solution of the linearized wave equation in a vacuum, they 
satisfy the Dirac gauge and are traceless (which is the linear approximation of unit determinant). 
The background is fiat, i.e., = ^ = 1, and /3* = 0. We display in Fig. [T]the radial profile of the 
component h'^^ at i = 6, at the equator and at the pole, with different values of the number of 
radial and angular points respectively. The analytical solution is recovered in the linear regime, 
as well as the velocity and the amplitude of the wave, and its decay with the radius. The 
absolute errors (L2 norm) in the numerical simulations are smaller than 8 x 10~^ for all the non 
zero components. The Sommerfeld condition at the outer boundary produces ingoing reflections 
that do not grow in time and have an amplitude much smaller than the amplitude of the initial 
one. We obtain second order of convergence for all the non zero components of the tensor h^^ , 
due to the influence of the inner boundary conditions imposed at r = 0, = and 9 = tt/2. 
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Figure 2. Scheme of the radial grid used in the code for rotating neutron star simulations. 



3.2. Equilibrium configuration of rotating neutron stars 

We consider an axisymmetric and uniformly rotating neutron star in equilibrium. The initial 
data have been obtained from LORENE [lOJ, a library of spectral methods. These models, of 
different rotation parameters, have non- vanishing /3^, N, ip and matter fields. All these variables 
are kept fixed during the evolution of the tensor h^^ . The initial model is a neutron star with a 
550 Hz rotation frequency, a 1.6 Mq baryon mass and a 12.86 km coordinate equatorial radius. 

Figure [2] shows a sketch of the grid, where the following domains have been defined: The 
matter domain (MD) contains the star. The propagation domain (PD) is the one where the 
waves are well-resolved. The region for extracting the gravitational radiation is placed at the 
outer boundary of the PD. The damping domain (DD) has a lower resolution in the radial 
direction. The wave travels far enough for imposing the Sommerfeld condition at the outer 
boundary before reaching it. Let tmb-pd be the radius which separates the MD and the PD, 
and rpD-DD be the radius which separates the PD and the DD. The star radius, r^, is close to 
th radius ^md-pd but still smaller in value. The logarithmic scale (in the radial direction) is a 
very useful tool for constructing this kind of grid with a reasonable number of points. We get 
the necessary accuracy in the MD and the PD, but the time stepping becomes more severe than 
the one required in an equally spaced radial grid. 

As initial data we have used stationary values for the tensor h'^^ and some differences from 
stationarity for A^^ . The aim at introducing these initial differences from stationarity is two- 
fold: on the one hand, the recovering of stationarity from a perturbed initial data, and, on the 
other hand, testing the outer boundary Sommerfeld condition imposed by the generation of the 
artificial wave. In this general background, those differences from stationarity introduced in 
the initial data, generate a perturbation that propagates to the outer boundary, leaving behind 
the stationary solution. In Fig. [3l we have plotted the absolute value of the component h^^ 
in terms of the radius for different times. We can see the perturbation traveling towards the 
outer boundary, where the wave leaves the numerical domain. The differences from stationarity 
are interpreted as an initial perturbation. The refiections close to the outer boundary come 
from the numerical implementation of the Sommerfeld condition. Second order of convergence 
is obtained for all the non zero components of the tensor h'^^ , as previously. 

This test helps to understand how the Sommerfeld condition at the outer boundary works. 
This condition is imposed to treat the outgoing waves which reach the outer boundary, as the 
one caused by the initial data. If the background is not fiat enough, the Sommerfeld condition 
can interpret the background as an outgoing wave. This behavior can be seen in Fig. H] for the 
simulation whose outer boundary is placed closer, at 3 x 10''. The Sommerfeld condition works 
properly for rout ^ 3 x lO'^cm (around 300 stellar radii). In our calculations we have used 80 
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Figure 3. Radial profile of \h''"^\. Dashed line corresponds to the stationary solution (as 
reference). Continuous lines correspond to the evolution of the initial data at 0.02, 0.05, 0.1, 1.0 
and 3.0 ms respectively. Vertical lines denote radius between domains. 
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Figure 4. Radial profile of \h'^'^\ at time t = 3 ms. The dashed line corresponds to the stationary 
solution (as reference). Continuous lines correspond to the numerical evolutions whose outer 
boundary are placed at 3 x 10'', 3 x 10*^ and 3 x 10^ cm. Vertical lines denote radii between 
domains. 



radial grid points in the MD, 80 points in the PD and around 50 in the DD. 

Once the outer boundary radius has been fixed, we can use the stationary initial data for 
both /i*-' and A^^ tensors. In Fig. [5] we show the absolute value of the component /i^'', in terms 
of the radius, for different times. No initial artificial wave is introduced. The noise comes from 
the outer boundary condition. Again, second order of convergence is obtained. 

3. 3. Perturbed equilibrium configuration of a rotating neutron star 

We consider a perturbed neutron star. The perturbation is proportional to sin(7r r/r^). Since 
the star is rotating, depends on 9, and the perturbation is not spherically symmetric; hence, 
gravitational radiation is expected to be generated. The hydrodynamic equations governing the 
matter evolution are solved. We keep fixed A^, /3* and ip, assuming that they do not change very 
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Figure 5. Radial profile of The dashed line corresponds to the stationary solution (as 

reference). Continuous lines correspond to the evolution of the initial data at 0.1, 1.0 and 3.0 
ms. Vertical lines denote radii between domains. 
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Figure 6. rhj^ versus the retarded time. The corresponding wave is extracted at r = 10 , 
r = 2.5 X 10^ and r = 5 x 10'' cm, at the equator, tpd-dd is placed at 3 x lO^cm. 



much during the evolution. This test will provide where to place rpD_DD; and, consequently, 
where to extract the gravitational radiation. 

In Fig.[6l an approximation of the real part of the Weyl scalar ^'4, (scaled with r) is plotted 
in terms of the retarded time, for different values of the extraction radius, at the equator. This 
approximation is accurate enough for the objectives of this test, tpd-dd is placed at 3 x lO^cm 
(around 30 stellar radii). The waves correspond to the physical gravitational waves coming from 
the evolution of the perturbed star. This evolution is governed by the coupled hydrodynamic 
and the Einstein equations (but keeping fixed A'^, /?* and ip). The different curves in Fig. [6] refer 
to different radii where the gravitational radiation has been extracted. From this figure, we can 
deduce that the speed of the waves is the light velocity and they decay as 1/r. We also see that 
the wave must be extracted far from the source but inside the PD. Regarding the resolution of 
the logarithmic grid, we have noticed that it is enough to cover a wavelength with five points at 
the end of the PD in order to obtain the wave accurately. 

Furthermore, it will be important to know where to place tpd-dd- In Fig- El fhj^ is plotted 
in terms of the retarded time in simulations with different tpd-dd- The wave is extracted close 
to rpD-DD in all the cases. From the figure, it is reasonable to place tpd-dd at 3 x 10''cm 
(around 30 stellar radii) or further. 
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Figure 7. r/i+ versus the retarded time, at the equator. Radii are expressed in cm. 
The corresponding wave is extracted at r = 10^ (?'pd-dd = 1-5 x 10''), r = 2.5 x 10^ 
(rpD-DD = 3 X 10"^) and r = 4.5 x 10^ (rpo-DD = 5 x 10^). 



4. Conclusions 

The first version of CoCoNuT's code [TT] was designed to evolve matter fields of a perfect 
fluid in the dynamical spacetime of the CFC approximation. The current version of CoCoNuT 
incorporates magnetic fields, i.e., it is a general-relativistic magneto-hydrodynamic code which 
evolves matter in the dynamical spacetime of CFC. The present contribution to this conference, 
together with previous works in [3, 4J, is a step towards the upgrading of the metric evolution 
of CoCoNuT's code by substituting the CFC approximation of the Einstein equations by the 
FCF of |2j. The hyperbolic sector to evolve h^^ has already been included into CoCoNuT's code. 
In this contribution we have shown that this sector works properly. We have presented the 
numerical evolution of Teukolsky waves in a flat background. We have analyzed the behavior 
of the code under the Sommerfeld boundary condition, and detected the outer radius, in the 
case of spacetimes generated by equilibrium configurations of rotating neutron stars. Finally, 
we have extracted the gravitational waveforms in the case of dynamical spacetimes generated 
by perturbed rotating neutron stars. 
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